###################################
## Refugees & Civil War          ##
## Functions Code                ##
###################################

## Load required packages:
require(arm)
require(lubridate)
require(stargazer)
require(lme4)
require(knitr)
require(reshape2)
require(estimatr)
require(tidyverse)
require(doMC)
require(texreg)
library(PanelMatch)
#library(fect)

require(ggplot2) #plotting packages
require(RColorBrewer) 
require(gridExtra)
require(stargazer)
require(xtable)
require(lattice)
require(grid)
require(patchwork)
require(ggcorrplot)

require(sp) #GIS packages
require(rgdal)
require(maps)
require(mapdata)
require(maptools)
require(ggmap)

## Not in 
'%!in%' <- function(x,y)!('%in%'(x,y))

## Shift lag and leads variables
shift <- function(x, shift_by){
  stopifnot(is.numeric(shift_by))
  stopifnot(is.numeric(x))
  
  if (length(shift_by)>1)
    return(sapply(shift_by,shift, x=x))
  
  out<-NULL
  abs_shift_by=abs(shift_by)
  if (shift_by > 0 )
    out<-c(tail(x,-abs_shift_by),rep(NA,abs_shift_by))
  else if (shift_by < 0 )
    out<-c(rep(NA,abs_shift_by), head(x,-abs_shift_by))
  else
    out<-x
  out
}

## Summary table
summaryfun <- function(x) as.data.frame(cbind(Mean = mean(x, na.rm = T), 
                                              Median = median(x, na.rm = T),
                                              SD = sd(x, na.rm = T),
                                              Min = min(x, na.rm = T),
                                              Max = max(x, na.rm = T),
                                              NAs = sum(is.na(x)) ))

## Transform BayesGLM class to GLM for stargazer
transbayesclass <- function(fit){
  class(fit) <- c("glm","lm")
  fit$call[1] <- quote(lm())
  return(fit)
}


## Predicted Probability functions ##

# for binary treatment (site) without interaction term
pred.bi <- function(fit){ #quasi-bayesian monte carlo 
  
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  reps <- matrix(NA, 2, 1000)
  
  for(i in 1:length(1:2)){
    if(i == 1){ #no site in my province 
      paste[,2] <- 0
    }else if(i == 2){ #site in my province
      paste[,2] <- 1
    }
    
    logistic <- function(x) exp(x)/(1+exp(x))
    yhat <- logistic(as.matrix(paste) %*% t(beta.sim))
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  means <- apply(reps, 1 , mean) 
  se <- apply(reps, 1 , sd)
  
  diff <- means[2] - means[1]
  diff.se <- sqrt(se[2]^2 + se[1]^2)
  
  means2 <- c(means[1], means[2], diff)
  
  se2 <- c(se[1], se[2], diff.se)
  
  ciu <- means2 + qnorm(0.975)*se2
  cil <- means2 - qnorm(0.975)*se2
  
  pred <- data.frame(Means = means2,
                     SE = se2,
                     CIU = ciu,
                     CIL = cil,
                     Est = c("No Refugee Presence",
                             "Refugee Presence",
                             "Effect of \nRefugee Presence"),
                     Plot = c(1:3))
  
  return(pred)
}

# for binary treatment (site) with interaction term
pred.bi.int <- function(fit){ #quasi-bayesian monte carlo 
  
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  reps <- matrix(NA, 4, 1000)
  
  for(i in 1:length(1:4)){
    if(i == 1){ #no site in my province or other province
      paste[,2] <- 0
      paste[,3] <- 0
    }else if(i == 2){ #site in my province, no site in other province
      paste[,2] <- 1
      paste[,3] <- 0
    }else if(i == 3){ #no site in my province, site in other province
      paste[,2] <- 0
      paste[,3] <- 1
    }else{ #site in my province, site in other province
      paste[,2] <- 1
      paste[,3] <- 1
    }
    paste[,4] <- paste[,2]*paste[,3]
    
    logistic <- function(x) exp(x)/(1+exp(x))
    yhat <- logistic(as.matrix(paste) %*% t(beta.sim))
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  means <- apply(reps, 1 , mean) 
  se <- apply(reps, 1 , sd)
  
  diff <- c(means[2] - means[1],
            means[4] - means[3])
  diff.se <- c(sqrt(se[2]^2 + se[1]^2),
               sqrt(se[4]^2 + se[3]^2))
  
  means2 <- c(means[1], means[2], diff[1],
              means[3], means[4], diff[2])
  
  se2 <- c(se[1], se[2], diff.se[1],
           se[3], se[4], diff.se[2])
  
  ciu <- means2 + qnorm(0.975)*se2
  cil <- means2 - qnorm(0.975)*se2
  
  pred <- data.frame(Means = means2,
                     SE = se2,
                     CIU = ciu,
                     CIL = cil,
                     Est = c("No Refugee \nPresence and \nNo Presence in \nOther Provinces",
                             "Refugee Presence but \nNo Presence in \nOther Provinces",
                             "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces",
                             "No Refugee Presence but \nPresence in \nOther Provinces",
                             "Refugee Presence and \nPresence in \nOther Provinces",
                             "Effect of Refugee \nPresence given \nPresence \nin Other Provinces"),
                     Plot = c(1:6))
  
  return(pred)
}

## Logistic regression with triple-interaction
pred.bi.tripint <- function(fit, hetvar, placebo = FALSE){
  
  ## Setup, draw coefficients
  logistic <- function(x) exp(x)/(1+exp(x))
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  ## Get indices
  if(placebo){
    maineff_rtb <- which(cols == "rtb.placebo")
  }else{
    maineff_rtb <- which(cols == "rtb") 
  }
  maineff_rtb.other <- which(cols == "rtb.other")
  maineff_hte <- which(cols == hetvar)
  
  if(placebo){
    inter_rtb_rtb.other <- which(cols == "rtb.placebo:rtb.other")
    inter_rtb_hte <- which(cols == paste0("rtb.placebo:", hetvar))
  }else{
    inter_rtb_rtb.other <- which(cols == "rtb:rtb.other")
    inter_rtb_hte <- which(cols == paste0("rtb:", hetvar))
  }
  inter_rtb.other_hte <- which(cols == paste0("rtb.other:", hetvar))
  
  if(placebo){
    tripint <- which(cols == paste0("rtb.placebo:rtb.other:", hetvar))
  }else{
    tripint <- which(cols == paste0("rtb:rtb.other:", hetvar)) 
  }
  
  ## ---------------------
  ## Concentrated, hte = 1
  ## ---------------------
  paste[,c(maineff_rtb.other, inter_rtb_rtb.other, 
           inter_rtb.other_hte, tripint)] <- 0
  
  paste[,c(maineff_rtb, maineff_hte, inter_rtb_hte)] <- 1
  yhat_rtb1 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  paste[,c(maineff_rtb, inter_rtb_hte)] <- 0
  yhat_rtb0 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  diff <- mean(yhat_rtb1) - mean(yhat_rtb0)
  diff.se <- sqrt(sd(yhat_rtb1)^2 + sd(yhat_rtb0)^2)
  
  df_1_0_1 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces", 
    hte = 1, hetvar = hetvar
  )
  
  ## ---------------------
  ## Concentrated, hte = 0
  ## ---------------------
  paste[,c(maineff_rtb.other, maineff_hte, 
           inter_rtb_rtb.other, inter_rtb_hte,
           inter_rtb.other_hte, tripint)] <- 0
  
  paste[,c(maineff_rtb)] <- 1
  yhat_rtb1 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  paste[,c(maineff_rtb)] <- 0
  yhat_rtb0 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  diff <- mean(yhat_rtb1) - mean(yhat_rtb0)
  diff.se <- sqrt(sd(yhat_rtb1)^2 + sd(yhat_rtb0)^2)
  
  df_1_0_0 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces", 
    hte = 0, hetvar = hetvar
  )
  
  ## ------------------  
  ## Dispersed, hte = 1
  ## ------------------
  paste[,c(maineff_rtb, maineff_rtb.other, maineff_hte,
           inter_rtb_rtb.other, inter_rtb_hte, inter_rtb.other_hte, 
           tripint)] <- 1
  yhat_rtb1 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  paste[,c(maineff_rtb, inter_rtb_rtb.other, 
           inter_rtb_hte, tripint)] <- 0
  yhat_rtb0 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  diff <- mean(yhat_rtb1) - mean(yhat_rtb0)
  diff.se <- sqrt(sd(yhat_rtb1)^2 + sd(yhat_rtb0)^2)
  
  df_1_1_1 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nPresence \nin Other Provinces", 
    hte = 1, hetvar = hetvar
  )
  
  ## ------------------
  ## Dispersed, hte = 0
  ## ------------------
  paste[,c(maineff_hte, inter_rtb_hte, 
           inter_rtb.other_hte, tripint)] <- 0
  
  paste[,c(maineff_rtb, maineff_rtb.other, inter_rtb_rtb.other)] <- 1
  yhat_rtb1 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  paste[,c(maineff_rtb, inter_rtb_rtb.other)] <- 0
  yhat_rtb0 <- apply(logistic(as.matrix(paste) %*% t(beta.sim)), 
                     2, mean, na.rm=T)
  
  diff <- mean(yhat_rtb1) - mean(yhat_rtb0)
  diff.se <- sqrt(sd(yhat_rtb1)^2 + sd(yhat_rtb0)^2)
  
  df_1_1_0 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nPresence \nin Other Provinces", 
    hte = 0, hetvar = hetvar
  )
  
  df_het_out <- rbind(df_1_0_1, df_1_0_0, df_1_1_1, df_1_1_0)
  df_het_out$Plot <- 1:nrow(df_het_out)
  return(df_het_out)
  
}

# for linear models with interaction term
pred.lm.int <- function(fit){ #quasi-bayesian monte carlo 
  
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  reps <- matrix(NA, 4, 1000)
  
  for(i in 1:length(1:4)){
    if(i == 1){ #no site in my province or other province
      paste[,2] <- 0
      paste[,3] <- 0
    }else if(i == 2){ #site in my province, no site in other province
      paste[,2] <- 1
      paste[,3] <- 0
    }else if(i == 3){ #no site in my province, site in other province
      paste[,2] <- 0
      paste[,3] <- 1
    }else{ #site in my province, site in other province
      paste[,2] <- 1
      paste[,3] <- 1
    }
    
    paste[,4] <- paste[,2]*paste[,3]
    
    yhat <- as.matrix(paste) %*% t(beta.sim)
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  means <- apply(reps, 1 , mean) 
  se <- apply(reps, 1 , sd)
  
  diff <- c(means[2] - means[1],
            means[4] - means[3])
  diff.se <- c(sqrt(se[2]^2 + se[1]^2),
               sqrt(se[4]^2 + se[3]^2))
  
  means2 <- c(means[1], means[2], diff[1],
              means[3], means[4], diff[2])
  
  se2 <- c(se[1], se[2], diff.se[1],
           se[3], se[4], diff.se[2])
  
  ciu <- means2 + qnorm(0.975)*se2
  cil <- means2 - qnorm(0.975)*se2
  
  pred <- data.frame(Means = means2,
                     SE = se2,
                     CIU = ciu,
                     CIL = cil,
                     Est = c("No Refugee \nPresence and \nNo Presence in \nOther Provinces",
                             "Refugee Presence but \nNo Presence in \nOther Provinces",
                             "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces",
                             "No Refugee Presence but \nPresence in \nOther Provinces",
                             "Refugee Presence and \nPresence in \nOther Provinces",
                             "Effect of Refugee \nPresence given \nPresence \nin Other Provinces"),
                     Plot = c(1:6))
  
  return(pred)
}

# for linear models robust with interaction term
pred.lm.int.r <- function(fit){ #lm robust
  
  ## When other = 0
  est_o0 <- coef(fit)[2]
  se_o0 <- sqrt(vcov(fit)[2,2])
  cil_o0 <- est_o0 - se_o0 * qnorm(.975)
  ciu_o0 <- est_o0 + se_o0 * qnorm(.975)
  
  ## When other = 1
  est_o1 <- coef(fit)[2] + coef(fit)[4]
  se_o1 <- sqrt(vcov(fit)[2,2] + vcov(fit)[4,4] + 2*vcov(fit)[2,4])
  cil_o1 <- est_o1 - se_o1 * qnorm(.975)
  ciu_o1 <- est_o1 + se_o1 * qnorm(.975)
  
  return(data.frame(Means = c(est_o0, est_o1), 
                    SE = c(se_o0, se_o1),
                    CIL = c(cil_o0, cil_o1),
                    CIU = c(ciu_o0, ciu_o1),
                    Est = c("Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces",
                             "Effect of Refugee \nPresence given \nPresence \nin Other Provinces"),
                    outcome = c(fit$outcome)))
  
}

pred.lm.tripint.r <- function(fit, hetvar, placebo = FALSE){
  
  ## Get indices
  cols <- names(coef(fit))
  coef_fit <- coef(fit)
  vcov_fit <- vcov(fit)
  if(placebo){
    maineff_rtb <- which(cols == "rtb.placebo")
  }else{
    maineff_rtb <- which(cols == "rtb") 
  }
  maineff_rtb.other <- which(cols == "rtb.other")
  maineff_hte <- which(cols == hetvar)
  
  if(placebo){
    inter_rtb_rtb.other <- which(cols == "rtb.placebo:rtb.other")
    inter_rtb_hte <- which(cols == paste0("rtb.placebo:", hetvar))
  }else{
    inter_rtb_rtb.other <- which(cols == "rtb:rtb.other")
    inter_rtb_hte <- which(cols == paste0("rtb:", hetvar))
  }
  inter_rtb.other_hte <- which(cols == paste0("rtb.other:", hetvar))
  
  if(placebo){
    tripint <- which(cols == paste0("rtb.placebo:rtb.other:", hetvar))
  }else{
    tripint <- which(cols == paste0("rtb:rtb.other:", hetvar)) 
  }
  
  ## ---------------------
  ## Concentrated, hte = 1
  ## ---------------------
  diff <- coef_fit[maineff_rtb] + coef_fit[inter_rtb_hte]
  diff.se <- sqrt(vcov_fit[maineff_rtb,maineff_rtb] + 
                    vcov_fit[inter_rtb_hte,inter_rtb_hte] + 
                    2*vcov_fit[maineff_rtb,inter_rtb_hte])
  df_1_0_1 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces", 
    hte = 1, hetvar = hetvar
  )
  
  ## ---------------------
  ## Concentrated, hte = 0
  ## ---------------------
  diff <- coef_fit[maineff_rtb]
  diff.se <- sqrt(vcov_fit[maineff_rtb,maineff_rtb])
  df_1_0_0 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nNo Presence \nin Other Provinces", 
    hte = 0, hetvar = hetvar
  )
  
  ## ------------------
  ## Dispersed, hte = 1
  ## ------------------
  diff <- coef_fit[maineff_rtb] + coef_fit[inter_rtb_hte] + 
    coef_fit[inter_rtb_rtb.other] + coef_fit[tripint]
  diff.se <- sqrt(
    vcov_fit[maineff_rtb,maineff_rtb] + 
      vcov_fit[inter_rtb_hte,inter_rtb_hte] + 
      vcov_fit[inter_rtb_rtb.other,inter_rtb_rtb.other] + 
      vcov_fit[tripint,tripint] + 
      2*vcov_fit[maineff_rtb,inter_rtb_hte] + 
      2*vcov_fit[maineff_rtb,inter_rtb_rtb.other] + 
      2*vcov_fit[maineff_rtb,tripint] + 
      2*vcov_fit[inter_rtb_hte,inter_rtb_rtb.other] + 
      2*vcov_fit[inter_rtb_hte,tripint] +
      2*vcov_fit[inter_rtb_rtb.other,tripint]
  )
  df_1_1_1 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nPresence \nin Other Provinces", 
    hte = 1, hetvar = hetvar
  )
  
  ## ------------------
  ## Dispersed, hte = 0
  ## ------------------
  diff <- coef_fit[maineff_rtb] + coef_fit[inter_rtb_rtb.other]
  diff.se <- sqrt(vcov_fit[maineff_rtb,maineff_rtb] + 
                    vcov_fit[inter_rtb_rtb.other,inter_rtb_rtb.other] + 
                    2*vcov_fit[maineff_rtb,inter_rtb_rtb.other])
  df_1_1_0 <- data.frame(
    Means = diff, SE = diff.se,
    CIU = diff + diff.se * qnorm(.975),
    CIL = diff - diff.se * qnorm(.975),
    Est = "Effect of Refugee \nPresence given \nPresence \nin Other Provinces", 
    hte = 0, hetvar = hetvar
  )
  
  df_het_out <- rbind(df_1_0_1, df_1_0_0, df_1_1_1, df_1_1_0)
  df_het_out$Plot <- 1:nrow(df_het_out)
  return(df_het_out)
  
}

# for two binary treatments (site) equivalent to categorical
pred.bi.cat <- function(fit){ #quasi-bayesian monte carlo 
  
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  reps <- matrix(NA, 3, 1000)
  
  for(i in 1:length(1:3)){
    if(i == 1){ #no old site refugees and no new site refugees
      paste[,2] <- 0
      paste[,3] <- 0
    }else if(i == 2){ #new site refugees, no old site refugees
      paste[,2] <- 1
      paste[,3] <- 0
    }else if(i == 3){ #no new site refugees, old site refugees
      paste[,2] <- 0
      paste[,3] <- 1
    }
    
    logistic <- function(x) exp(x)/(1+exp(x))
    yhat <- logistic(as.matrix(paste) %*% t(beta.sim))
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  means <- apply(reps, 1 , mean) 
  se <- apply(reps, 1 , sd)
  
  diff <- c(means[2] - means[1],
            means[3] - means[1])
  diff.se <- c(sqrt(se[2]^2 + se[1]^2),
               sqrt(se[3]^2 + se[1]^2))
  
  means2 <- c(diff[1], diff[2])
  
  se2 <- c(diff.se[1], diff.se[2])
  
  ciu <- means2 + qnorm(0.975)*se2
  cil <- means2 - qnorm(0.975)*se2
  
  pred <- data.frame(Means = means2,
                     SE = se2,
                     CIU = ciu,
                     CIL = cil,
                     Est = c("Effect of New Refugee Presence",
                             "Effect of Protracted Refugee Presence"),
                     Plot = c(1:2))
  
  return(pred)
}

# for linear models with two binary treatments (site) equivalent to categorical
pred.lm.cat <- function(fit){ #quasi-bayesian monte carlo 
  
  beta.sim <- as.data.frame(mvrnorm(1000, mu = coef(fit), Sigma = vcov(fit)))
  paste <- model.matrix(fit)
  cols <- colnames(paste)
  beta.sim <- beta.sim[, c(cols)]
  
  reps <- matrix(NA, 3, 1000)
  
  for(i in 1:length(1:3)){
    if(i == 1){ #no old site refugees and no new site refugees
      paste[,2] <- 0
      paste[,3] <- 0
    }else if(i == 2){ #new site refugees, no old site refugees
      paste[,2] <- 1
      paste[,3] <- 0
    }else if(i == 3){ #no new site refugees, old site refugees
      paste[,2] <- 0
      paste[,3] <- 1
    }
    
    yhat <- as.matrix(paste) %*% t(beta.sim)
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  
  means <- apply(reps, 1 , mean) 
  se <- apply(reps, 1 , sd)
  
  diff <- c(means[2] - means[1],
            means[3] - means[1])
  diff.se <- c(sqrt(se[2]^2 + se[1]^2),
               sqrt(se[3]^2 + se[1]^2))
  
  means2 <- c(diff[1], diff[2])
  
  se2 <- c(diff.se[1], diff.se[2])
  
  ciu <- means2 + qnorm(0.975)*se2
  cil <- means2 - qnorm(0.975)*se2
  
  pred <- data.frame(Means = means2,
                     SE = se2,
                     CIU = ciu,
                     CIL = cil,
                     Est = c("Effect of New Refugee Presence",
                             "Effect of Protracted Refugee Presence"),
                     Plot = c(1:2))
  
  return(pred)
}

# for linear models robust with two binary treatments (site) equivalent to categorical
pred.lm.cat.r <- function(fit){ #lm robust
  
  ## When new site
  est_o0 <- coef(fit)[2]
  se_o0 <- sqrt(vcov(fit)[2,2])
  cil_o0 <- est_o0 - se_o0 * qnorm(.975)
  ciu_o0 <- est_o0 + se_o0 * qnorm(.975)
  
  ## When old site 
  est_o1 <- coef(fit)[3]
  se_o1 <- sqrt(vcov(fit)[3,3])
  cil_o1 <- est_o1 - se_o1 * qnorm(.975)
  ciu_o1 <- est_o1 + se_o1 * qnorm(.975)
  
  return(data.frame(term = names(fit$coefficients[c(2,3)]),
                    estimate = c(est_o0, est_o1), 
                    std.error = c(se_o0, se_o1),
                    conf.low = c(cil_o0, cil_o1),
                    conf.high = c(ciu_o0, ciu_o1),
                    Est = c("Effect of New Refugee Presence",
                            "Effect of Protracted Refugee Presence"),
                    outcome = c(fit$outcome)))
  
}

